%This program localizes single molecules and determines their positions
%by fitting a gaussian to their intensity profile.  A position dependent
%background subtraction is used.

function einzelreader2(handles, h)
    global xpix ypix wbox;

    version = '6/10/08'; %record version (date) of gui being used

    rbox=str2double(get(handles.rbox_edit,'String'));
    wbox=2*rbox+1;
    [xpix,ypix] = meshgrid(-rbox:rbox,-rbox:rbox);
    
    q       = str2double(get(handles.cam_pix_edit,'String')); 
    wvlnth  = str2double(get(handles.wvlnth,'String'))/1000; %convert wavelength from nm to um       
    NA      = str2double(get(handles.NA,'String'));
    psf_scale = str2double(get(handles.psf_scale_edit,'String'));
    
    psf_w0 = psf_scale*0.55*wvlnth/NA/1.17; % 1/e2 radius of PSF in um, use 1/e2 = FWHM/1.17 from Pawley
                                      % with scale factor (20% for "real objective" due to
                                      % measured PSF from Hess and Webb
                                      % 2002)
    psf_std=psf_w0/2; %standard deviation of psf
    psf_w02=(psf_w0/q)*(psf_w0/q); %square of 1/e^2 radius in pixels
    
    yfit_psf=exp(-2*((xpix).*(xpix)+(ypix).*(ypix))/psf_w02);
    npix=sum(sum(yfit_psf));     % area of molecule in square pixels
    
    total_molecules=0;
    xcm_all=zeros(1,1);
    ycm_all=zeros(1,1);
    framenum_all=zeros(1,1);
    n_fail_a0=0;
    n_fail_outbox=0;

    imagefile = get(handles.im_file_edit, 'String');
    outpath = get(handles.out_dir_edit, 'String');
    if outpath(length(outpath))~='\'
        outpath = [outpath,'\'];
    end
    preview = get(handles.show_preview,'Value');

    fs0 = str2double(get(handles.fs,'String')); %Setup for preview frame skipping
    fs1 = fs0;
    
    if get(handles.rolling_ball_cb,'Value') %if using rolling ball subtraction
        %Obtain info on image type, etc
        file_info = get_file_info2(imagefile, handles);
        
        n_end   = file_info.stop;
        n_start = file_info.start;
        
        if get(handles.custom_roi,'Value')
            x_offset   = str2double(get(handles.x_off_edit,'String'));
            y_offset   = str2double(get(handles.y_off_edit,'String'));
            x_size     = str2double(get(handles.x_size_edit,'String'));
            y_size     = str2double(get(handles.y_size_edit,'String'));
        else
            x_offset = 1;
            y_offset = 1;
            x_size   = file_info.width;
            y_size   = file_info.height;
        end
        
        rball=str2double(get(handles.rb_radius_edit,'String')); %radius of rolling ball
        se = strel('ball',rball,rball,0); %structural element, i.e. rolling ball
        
        FWHM=str2double(get(handles.FWHM_edit,'String'));; %FWHM of gaussian smoothing in pixels
        rk=(FWHM)/sqrt(2*log(2)); %1/e^2 smoothing radius in pixels
        kw=20; %kernal width of smoothing function
        [X,Y]=meshgrid(-kw/2:kw/2,-kw/2:kw/2);
        kd=sqrt(X.*X+Y.*Y);
        gs=exp(-2*kd.*kd/(rk*rk));
        gs=gs/sum(sum(gs)); %smoothing function normalized to have area = 1
        
    else %else use integrated widefield sum
        %Load wf sum & mean background mat file
        meanbkg_mat_file = getappdata(0,'last_mbkg_save'); %strcat(dir,'analysis\',date,'\',base_name,'sum_wf_meanbkg.mat');
        load(meanbkg_mat_file,'meanbkg_all','image_sum','x_size','y_size','x_offset','y_offset','field0','n_field','file_info');%,'zero_level'); %mat file from sum_wf_meanbkg.m
        
        n_end   = file_info.stop;
        n_start = file_info.start;
        %give warning if starting/ending frames are not compatible with wf sum & mean background mat file 
        if(n_start<field0)
            warndlg(sprintf('Warning: Starting frame is less than corresponding frame used in Compute WF & Mean Background, need to change starting frame or re-run Compute WF & Mean Background'));
            return  
        end
        if(n_end>n_field)
            warndlg(sprintf('Warning: Ending frame is greater than corresponding frame used in Compute WF & Mean Background, need to change ending frame or re-run Compute WF & Mean Background'));
            return  
        end
        
        bkg_pc = str2double(get(handles.bkg_percent,'String'))/100;
        zero_level = str2double(get(handles.zero_lvl,'String'));
        
        %normalize so that ave pixel value = 1
        image_sum_mean=mean(mean(image_sum(y_offset:y_offset+y_size-1,x_offset:x_offset+x_size-1)));
        image_sum_norm = (image_sum(y_offset:y_offset+y_size-1,x_offset:x_offset+x_size-1))/image_sum_mean;
    end
    
    %Load variables from the GUI edit boxes
    %--------------------------------------
    if get(handles.pixel_photon,'Value')
        pix_to_pho = str2double(get(handles.pix_to_pho,'String'));
    else
        pix_to_pho = 1;
    end

    iprod_thresh  = str2double(get(handles.iprod_thresh_edit,'String'))*pix_to_pho; %minimum threshold for a "bright object"    
    box_overlap_factor = str2double(get(handles.box_overlap_edit,'String')); %if center to center closer, don't include either
    w_mask = round(rbox*box_overlap_factor); %width of masking when a "high" pixel is identified
    
    %Initialize temporary runtime file for passing variables
    %-------------
    for x=1:99
        temp_runtime = [outpath,'temp_runtime_',num2str(x,'%2.2d'),'.mat'];
        if ~exist(temp_runtime,'file')
            break
        end
    end
    setappdata(0,'temp_runtime',temp_runtime);
    %-------------

    %Initialize waitbar and pause component
    w = waitbarxmod(0,'Executing "einzelreader.m" ...','CreateCancelBtn','delete(gcf)');
    set(w,'Name','Progress Bar');
    uicontrol('Style','pushbutton','Parent',w,'String','Pause','Position',[210,10,60,23], ...
              'UserData',1,'Callback',@pause_gui);
    drawnow;   %Draw the extra button immediately    
    pause(.1); %Pause to ensure window completes drawing
    %keepontop('Progress Bar');

    wb_norm = n_end-n_start;
    if wb_norm == 0 %Avoid divide by zero errors
        wb_norm = 1;
    end

    if h~=0                             %If preview window is active, store the handle of the progress bar
        setappdata(h,'PrBar_Handle',w); %in the appdata of the figure
    end

    for fileloop=n_start:n_end
        if strcmp(file_info.type,'singlepage')
            index = num2str(fileloop,file_info.prec);
%             infile = [file_info.path,'\',file_info.part_name,index,file_info.ext];
            infile = [file_info.path,'/',file_info.part_name,index,file_info.ext];
            i1=double(imread(infile));
        elseif strcmp(file_info.type,'multipage')
            i1=double(imread(imagefile,fileloop));
        end

        compl = (fileloop-n_start)/wb_norm;
        drawnow;
        waitbarxmod(compl,w,'Executing "einzelreader2.m" ...'); %Rename the waitbar
        
        if get(handles.rolling_ball_cb,'Value') %if using rolling ball subtraction
            i1_gs = uint16(conv2(i1,gs,'same')); %smoothed original
            rb_im = imopen(i1_gs,se);
            bkg   = zeros(y_size,x_size);
            bkg   = double(rb_im(y_offset:y_offset+y_size-1,x_offset:x_offset+x_size-1));  
        else
            bkg    = zeros(y_size,x_size);
            bkg    = image_sum_norm*bkg_pc*meanbkg_all(fileloop)+zero_level;
        end

        iprod = zeros(y_size,x_size);
        iprod = i1(y_offset:y_offset+y_size-1,x_offset:x_offset+x_size-1)-bkg;
        iprod=iprod.*(iprod>0); % set any negative pixel values to zero

        iprod_copy=iprod;
        high_pixel_mask = zeros(y_size,x_size);
        n_boxes=0;
        boxes_xy=zeros(10000,3);
        pix_val=iprod_thresh;

        while(pix_val >= iprod_thresh) %continue until reach minimum threshold
            pix_val=0;
            for i=rbox+1:y_size-rbox-1
                for j=rbox+1:x_size-rbox-1
                    if(high_pixel_mask(i,j)==0 && iprod_copy(i,j)>pix_val)
                        pix_val=iprod_copy(i,j); %find brightest pixel
                        high_pixel_y=i;
                        high_pixel_x=j;
                    end
                end
            end
            if(pix_val < iprod_thresh)
                break
            end
        
            x0_box=high_pixel_x-rbox;
            y0_box=high_pixel_y-rbox;
            x1_box=high_pixel_x+rbox;
            y1_box=high_pixel_y+rbox;
            
            x0_mask=high_pixel_x-w_mask;
            if x0_mask < 1
                x0_mask = 1;
            end
            x1_mask=high_pixel_x+w_mask;
            if x1_mask > x_size
                x1_mask = x_size;
            end
            y0_mask=high_pixel_y-w_mask;
            if y0_mask < 1
                y0_mask = 1;
            end
            y1_mask=high_pixel_y+w_mask;
            if y1_mask > y_size
                y1_mask = y_size;
            end
        
            high_pixel_mask(y0_mask:y1_mask,x0_mask:x1_mask)=1;

            grab=iprod_copy(y0_box:y1_box,x0_box:x1_box);
            grab_stdev=sqrt(grab);
            minbkg=min(min(grab)); 
            grab=grab-minbkg;
            grab_sum=sum(sum(grab));
        
            %calculate C.O.M. for intial guess in fit
            xm_sum=0;
            ym_sum=0;
            m_sum=0;
            for i=x0_box:x1_box
                for j=y0_box:y1_box
                    xind=floor(i);
                    yind=floor(j);
                    intens=iprod_copy(yind,xind);
                    xm_sum=xm_sum+xind*intens;
                    ym_sum=ym_sum+yind*intens;
                    m_sum=m_sum+intens;
                end
            end
            x_cm=xm_sum/m_sum;
            y_cm=ym_sum/m_sum;

            xc_box=(x0_box+x1_box)*0.5;
            yc_box=(y0_box+y1_box)*0.5;

            xguess=x_cm-xc_box;
            yguess=y_cm-yc_box;

            for i=1:wbox
                for j=1:wbox
                    k=(i-1)*wbox+j;
                    xymerge(k)=0;
                    zmerge(k)=grab(i,j);
                end
            end

            beta=[xguess,yguess,50,psf_w0/q];
            [betafit,resid,J,COVB,mse] = nlinfit(xymerge,zmerge,@gaussian_merge2,beta);
            ci = nlparci(betafit,resid,'covar',COVB); %calculate error estimates on parameters
            ci_err=(ci(:,2)-ci(:,1))/2;

            zfitmerge=gaussian_merge2(betafit,xymerge);

            for i=1:wbox
                for j=1:wbox
                    k=(i-1)*wbox+j;
                    zfit(i,j)=zfitmerge(k);
                end
            end

            yf=betafit(2)+yc_box;
            xf=betafit(1)+xc_box;
            a0=betafit(3);
            r0=abs(betafit(4));
        
            failed=0; %flag for failed localization
            if(a0 < 0)
                n_fail_a0=n_fail_a0+1;
                failed=1;
            end
            if(xf > x1_box || xf < x0_box || yf > y1_box || yf < y0_box)
                n_fail_outbox=n_fail_outbox+1;
                failed=1;
            end
            
            if(failed==0) %assign if fit criteria is satisfied
                total_molecules=total_molecules+1;
                xcm_all(total_molecules)=x_cm;
                ycm_all(total_molecules)=y_cm;
                xf_all(total_molecules)=xf;
                yf_all(total_molecules)=yf;
                a0_all(total_molecules)=a0;
                r0_all(total_molecules)=r0;
                framenum_all(total_molecules)=fileloop;
                pix_val_all(total_molecules)=pix_val;
                xf_err_all(total_molecules)=ci_err(1);
                yf_err_all(total_molecules)=ci_err(2);
                a0_err_all(total_molecules)=ci_err(3);
                r0_err_all(total_molecules)=ci_err(4);
                grab_sum_all(total_molecules)=grab_sum;

                n_boxes=n_boxes+1;
                boxes_xy(n_boxes,1)=high_pixel_x;
                boxes_xy(n_boxes,2)=high_pixel_y;
                boxes_xy(n_boxes,3)=1;
            end

            nmol_all(fileloop)=n_boxes;

            drawnow   %Flush event queue to keep cancel/pause button responsive on every pass.
                      %This is the slowest function by far and drawnow barely changes its exec time

            waitbarxmod(compl,w); %check waitbar 
        end
        
        if preview  %If preview active, update images on screen
            set(0,'CurrentFigure',h);
            if ~exist('plot2','var')
                position2 = [.01,.37,.25,.25];
                plot2 = axes('Position', position2);
                setappdata(plot2,'pos',position2);
            else
                set(h,'CurrentAxes',plot2);
            end
            i1display = i1(y_offset:y_offset+y_size-1,x_offset:x_offset+x_size-1);
            imagesc(i1display,'HitTest','off');
            axis image;
            colormap gray;
            title(['Last Frame Processed: ', num2str(fileloop)])
            set(plot2, 'ButtonDownFcn',{@focus_swap,guidata(gcbo)});
            clear i1display;

            if ~exist('plot1','var')
                position1 = [.01,.05,.25,.25];
                plot1 = axes('Position', position1,'Color','k','YTickLabel','','XTickLabel','');
                axis image;
                set(plot1,'Position',position1);
                setappdata(plot1,'pos', position1);
            end

            if get(handles.rolling_ball_cb,'Value')
                if ~exist('plot3','var')
                    position3 = [.01,.70,.25,.25];
                    plot3 = axes('Position', position3);
                    setappdata(plot3,'pos',position3);
                else
                    set(h,'CurrentAxes',plot3);
                end
                imagesc(rb_im,'HitTest','off')
                colormap gray, axis image;
                title(['Rolling Ball Profile: ',num2str(fileloop)])
                set(plot3,'ButtonDownFcn',{@focus_swap,guidata(gcbo)});
            end

            if ~exist('plot4','var')
                if x_size < y_size
                    positionf = [.18,.065,.87,.87];
                else
                    positionf = [.32,.3,.6,.65];
                end
                plot4 = axes('Position', positionf);
                setappdata(0,'focused',plot4);
            else
                set(h,'CurrentAxes',plot4);
            end

            if get(handles.mnl_scale_check,'Value')
                clim = str2double(get(handles.mnl_scale_edit,'String'));
                imagesc(iprod,'HitTest','off',[0 clim]);
            else
                imagesc(iprod,'HitTest','off');
            end

            axis image;
            colormap gray;
            title(['Einzel Reader: Frame: ',num2str(fileloop)]);
            set(plot4, 'ButtonDownFcn',{@focus_swap,guidata(gcbo)});

            if (get(handles.show_colorbar,'Value') == 1)
                focused = getappdata(0, 'focused');
                set(h,'CurrentAxes', focused);
                pos = get(focused, 'Position');
                colorbar('Position',[.947, pos(2), .015, pos(4)],'DrawMode','fast')
            end

            set(h,'CurrentAxes', plot4);
            hold on
            draw_boxes(n_boxes,boxes_xy,rbox);
            drawnow;
            hold off
            
            %If preview active, compute wf & merge on first, last, and nth frames (always)
            if fileloop == n_end
                setappdata(0,'Save',1); %Trigger image save on last frame
                save(temp_runtime);     %Save variables to a temporary file for function passing
                fpalm_render_einzelreader(handles,h);
            elseif fileloop == n_start+fs1-1
                fs1=fs1+fs0;
                save(temp_runtime);
                fpalm_render_einzelreader(handles,h);
            elseif fileloop == n_start
                save(temp_runtime);
                fpalm_render_einzelreader(handles,h);
            end
        end 
        
        drawnow;              %Draw any boxes and update einzelreader image (flush queue and check waitbar if preview not active)
        waitbarxmod(compl,w); %Update the waitbar
        
        if get(handles.frame_delay_check,'Value')
            delay = str2double(get(handles.frame_delay_edit,'String'))/1000;
            if delay ~= inf
                pause(delay);
            end
        end
    end

    if ~preview %If preview is off, continue to wf & merge at end and save the result
        plot1 = 0;
        save(temp_runtime);
        setappdata(0,'Save',1);
        fpalm_render_einzelreader(handles,h);
    else        %Otherwise delete progress bar and handle stored in appdata
        delete(w);
        rmappdata(h,'PrBar_Handle');
    end
    
    if get(handles.oc_wf,'Value') %save widefield image if selected
        if (file_info.part_name(length(file_info.part_name)) == '_') %Avoid saving with repeated seperators                
            wf_file = [outpath,'\',file_info.part_name,'wf.tif'];
        else
            wf_file = [outpath,'\',file_info.part_name,'_wf.tif'];
        end
        imwrite(uint8((image_sum/max(image_sum(:)))*255),wf_file,'Compression','none');
    end

%Save Variables
%----------------
    if (file_info.part_name(length(file_info.part_name)) == '_') %Avoid saving with repeated seperators
        out_file = [outpath,file_info.part_name,num2str(n_start),'-',num2str(n_end),'_t',num2str(iprod_thresh),'.mat'];
    else
        out_file = [outpath,file_info.part_name,'_',num2str(n_start),'-',num2str(n_end),'_t',num2str(iprod_thresh),'.mat'];
    end

    if ~exist(outpath,'dir') %Make sure the directory exists
        answer = questdlg('Output directory could not be accessed. Choose a new save location?','Warning!','Yes','Discard','Yes');

        if strcmp(answer,'Yes')
            outpath = uigetdir();
            if ~outpath
                return
            end
            if outpath(length(outpath))~='\'
                outpath = [outpath,'\'];
            end

            if (file_info.part_name(length(file_info.part_name)) == '_')
                out_file = [outpath,file_info.part_name,num2str(n_start),'-',num2str(n_end),'_t',num2str(iprod_thresh),'.mat'];
            else
                out_file = [outpath,file_info.part_name,'_',num2str(n_start),'-',num2str(n_end),'_t',num2str(iprod_thresh),'.mat'];
            end
        else
            return
        end
    end

    prec_edit     = get(handles.prec_edit,'String');                %<- These to be saved with the
    bkgn          = str2double(get(handles.bkgn_noise,'String'));   %<- mat file for loading purposes    
    ppp           = str2double(get(handles.pix_to_pho,'String'));   %<-
    jshift        = str2double(get(handles.jshift,'String'));       %<-
    ishift        = str2double(get(handles.ishift,'String'));       %<-
    exf           = str2double(get(handles.exf_edit,'String'));     %<-   expansion factor
    
    a0_phot=a0_all/ppp;    % peak amplitude of each molecule in photons
    N=npix*a0_phot;        % number of photons for each molecule

    %localization precision in um
    lp2=((psf_std^2)+(q^2)/12)*1./N+8*pi*(psf_std^4)*(bkgn^2)/(q^2)*1./(N.*N);
    lp=sqrt(lp2);
    
    list = who; %Use regular expressions to ensure handles are not saved
    match = regexp(list,'^plot|handles|\<h\>|\<w\>');
    for i=length(match):-1:1
        if match{i}
            list(i)=[];
        end
    end
    save(out_file,list{:}); %Save variables used in einzelreader.m
    setappdata(0,'last_einzel_save',out_file);  %Save location of einzelreader mat file for other functions to use
%----------------
    %Clean up temp files and appdata
    if exist(temp_runtime,'file')
        delete(temp_runtime)
    end
    if isappdata(0,'temp_runtime')
        rmappdata(0,'temp_runtime')
    end